Quarterly  Journal  of  the  Royal  Meteorological  Society 


Q.  J.  R.  Meteorol.  Soc.  00:  1  16  (2013) 


Grid  effect  on  spherical  shallow  water  jets  using 
continuous  and  discontinuous  Galerkin  methods 


S.  Marias*,  M.  A.  Kopera,  F.  X.  Giraldo 


Department  of  Applied  Mathematics,  Naval  Postgraduate  School 

*  Correspondence  to:  833  Dyer  Road,  Bldg.  232,  Spanagel  249A,  Monterey,  CA  93943  U.S.A. 

smarrasl@nps.edu;  fxgirald@nps.edu 


Element-based  Galerkin  methods  are  not  constrained  by  the  logical 
structure  of  the  computational  grid.  In  this  paper,  the  behavior  of  the 
continuous  and  discontinuous  Galerkin  solutions  of  the  shallow  water 
equations  on  the  sphere  is  analyzed  for  three  different  grids  of  ubiquitous  use 
in  atmospheric  modeling:  (A)  hexahedron  (or  cubed  sphere),  (B)  reduced 
latitude-longitude,  and  (G)  icosahedron.  Conforming  and  nonconforming 
mesh  configurations  are  adopted.  The  nonconforming  grids  are  based  on 
a  quad-tree  structure.  The  analysis  is  performed  on  the  mid-latitude  zonal 
flow  problem  suggested  by  Galewsky  et  al.  [An  initial-value  problem  for  testing 
numerieal  models  of  the  global  shallow-water  equation,  Tellus  2004;  56A:429-440]. 
This  test  is  sufficiently  simple  in  its  implementation,  yet  sufficiently  complex 
to  capture  some  important  modes  that  arise  on  the  global  scales  of  real 
atmospheric  dynamics.  Because  the  inviscid  solution  on  certain  grids  shows 
a  high  sensitivity  to  the  resolution,  the  viscous  counterpart  of  the  governing 
equations  are  also  solved  and  the  results  are  compared.  The  study  shows 
that  not  only  the  resolution,  but  the  grid  alignment  with  respect  to  the  flow 
is  important  in  order  to  obtain  an  accurate  solution.  This  is  especially  true  if 
the  governing  equations  are  not  regularized  by  the  addition  of  a  sufficiently 
large,  fully  artificial,  diffusion  mechanism;  we  give  a  brief  (not  exhaustive) 
analysis  of  the  effect  of  viscosity  on  the  solution  quality.  The  flexibility  and 
accuracy  of  element-based  Galerkin  methods  on  general  spherical  geometries 
is  illustrated;  we  particularly  emphasize  the  excellent  properties  of  the 
reduced  Lat-Lon  configuration  in  comparison  with  the  cubed-sphere  on  the 
one  hand,  and  with  the  more  regular  and  uniform  icosahedral  grid  on  the 
other. 
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As  computational  power  increases,  meteorologists  demand 
greater  resolution  from  global  atmospheric  models  utilizing 
spherical  domains.  Driven  by  the  need  to  select  the 
most  proper  computational  grid  on  the  sphere  for  the 
Nonhydrostatic  Unified  Model  of  the  Atmosphere  (NUMA) 
(Kelly  and  Giraldo  2012;  Giraldo  et  al.  2013),  in  this  paper 
we  analyze  the  continuous  Galerkin  (or  spectral  elements 


1.  Introduction 


from  now  on)  and  discontinuous  Galerkin  (DG)  solutions 
of  the  shallow  water  equations  (SWE)  on  a  set  of  common 
grids  used  in  global  circulation  models  (GCMs).  Specihcally, 
we  study  1)  the  cubed-sphere  (Hex)  (Sadourny  1972), 
2)  the  icosahedron  (Ico)  (Sadourny  et  al.  1968),  and  3) 
a  reduced  latitude-longitude  (Lat-Lon)  geometry  (Phillips 
1957).  Because  SWE  capture  many  of  the  essential  features 
of  GCMs  while  eliminating  the  vertical  structure  of  the 
atmosphere,  the  results  from  this  study  will  be  directly 
applicable  to  the  fully  three-dimensional  model  NUMA. 


t  Please  ensure  that  you  use  the  most  up  to  date  class  file,  available  from 
the  QJRMS  Home  Page  at 
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The  logically  structured  Hex  and  Lat-Lon  grids  can 
be  advantageous  for  a  hnite  difference  based  solver,  and 
are  generally  usable  by  grid  point  methods  such  as  finite 
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and  spectral  element  (FE,  SE),  discontinuous  Galerkin 
(DG),  or  finite  volume  (FV)  methods.  Furthermore,  Lat- 
Lon  is  the  standard  discretization  for  GGMs  based  on 
spectral  transform  methods  due  to  the  fast  Fourier  transform 
operation  along  the  longitude  direction  (Hogan  et  al.  1991). 
Icosahedral  grids  or,  more  generally,  unstructured  tri-  and 
quad-based  tessellations*  of  the  sphere  are  geometrically 
flexible,  but  not  all  numerical  methods  are  able  to  handle 
them.  (A  review  of  different  grid  generation  techniques 
in  atmospheric  modeling  can  be  found  in  the  book  by 
Behrens  (2006).)  Recently,  different  grids  to  be  used 
with  finite  volume  discretization  techniques  have  been 
analyzed  in,  e.g.,  Qaddouri  et  al.  (2012);  Weller  et  al.  (2012) 
or  Peixoto  and  Barros  (2013),  with  a  review  appearing 
in  Staniforth  and  Thuburn  (2012).  In  the  context  of 
Galerkin  methods,  however,  not  much  analysis  has  been 
conducted  with  respect  to  the  solution  dependence  on  these 
computational  grids.  Rather,  given  a  specihc  grid,  we  are 
likely  to  find  a  model  that  is  developed  around  it  and  is 
optimized  to  minimize  the  error  in  a  specihc  conhguration. 
In  the  case  of  the  solution  of  shallow  water  problems, 
examples  of  this  approach  are  the  works  by  Taylor  et  al. 
(1997)  (SE  on  cubed-sphere),  Lanser  et  al.  (2000)  (FV  on  a 
reduced  Lat-Lon  geometry),  Giraldo  (2001)  (SE  on  the  quad- 
based  icosahedron),  Giraldo  et  al.  (2002)  (DG  on  a  triangle- 
based  icosahedron),  or  Nair  et  al.  (2007)  (DG  on  the  cubed 
sphere).  One  step  forward  in  grid  comparison  was  made 
by  Giraldo  and  Warburton  (2005)  who  compared  different 
triangle-based  unstructured  grids  amongst  each  other  and, 
subsequently,  against  the  cubed  sphere.  They  concluded  that 
the  accuracy  of  the  solution  is  clearly  a  function  of  the 
combination  of  the  grid  and  the  problem  that  is  being  solved, 
and  it  does  not  depend  on  either  the  numerical  method  or 
the  grid  alone.  A  comparison  between  dynamically  adaptive 
and  uniform  meshes  was  performed  by  Muller  et  al.  (2013) 
for  triangle  based  discontinuous  Galerkin  methods. 

In  2008,  St-Cyr  et  al.  compared  two  shallow  water  solvers 
based  on  spectral  elements  on  the  one  hand,  and  hnite 
volumes  on  the  other  (see  St-Cyr  et  al.  (2008)).  The  two 
models  were  designed  to  execute  on  different  grids:  a  cubed- 
sphere  (SE  solver)  and  on  a  Lat-Lon  grid  (FV  model).  In 
their  study,  they  assert  that  the  spectral  element  method 
is  unable  to  capture  the  solution  at  certain  resolutions. 
In  this  paper  we  show  that  the  inaccuracies  that  they 
encountered  are  not  linked  to  the  numerical  method  per  se, 
but  rather,  to  the  unfortunate  combination  of  the  resolution 
and  alignment  of  the  mesh  with  respect  to  the  dynamics  that 
was  being  resolved  in  that  specific  study.  We  will  confirm 
their  conclusion  that  a  coarse  resolution  is  responsible  for 
a  wavenumber-4  signal  that  destroys  the  solution  on  the 
hexahedral  grid;  nevertheless,  the  simple  change  of  grid  will 
disproof  the  responsibility  of  SE  and,  rather,  show  that 
the  instabilities  are  merely  an  undesired  effect  of  the  grid 
conhguration. 

This  study  is  performed  by  solving  the  two  nonlinear 
zonal  hows  proposed  by  Galewsky  et  al.  (2004).  At  low 
resolution,  the  hows  on  the  cubed-sphere  and  icosahedron 
completely  diverge  from  the  true  solutions  in  the  case  of  (i) 
the  equilibrium,  and  (ii)  the  barotropic  instability  tests.  As 
the  resolution  is  increased,  the  effect  of  the  irregular  grid 
geometry  is  partially  camouhaged  and  eventually  disappears. 
In  search  of  the  coarsest  allowable  resolution  to  be  used  while 
still  preserving  accuracy,  we  added  an  artihcial  diffusion 


*The  keywords  grid,  mesh,  and  tessellation  will  be  used  interchangeably 
throughout  the  paper. 


with  constant  coefficient  i/ =  1  x  10®  m^  s“^  as  suggested 
by  Galewsky  et  al.  (2004).  The  positive  effect  of  diffusion 
is  clear  for  the  low  resolution  results  on  the  cubed-sphere 
and  icosahedron;  on  the  contrary,  we  will  show  how  this 
is  not  completely  true  for  the  solution  on  the  reduced 
Lat-Lon  grid.  It  must  be  pointed  ont  that  the  diffusion 
correction  that  was  adopted  is  fully  artificial  and  has  no 
physical  meaning.  The  sole  reason  for  its  use  in  this  study 
is  to  be  able  to  analyze  the  solution  when  the  grid  is 
not  of  sufficient  resolution  to  properly  resolve  the  flow. 
Partially  based  on  this  study,  but  mostly  on  the  analysis 
made  by  practitioners  in  computational  fluid  dynamics  for 
more  than  three  decades,  we  strongly  believe  that  simple 
Laplace  operators  (of  any  order)  for  stabilization  purposes 
should  be  avoided  for  the  sake  of  accuracy.  If  viscosity  is 
indeed  necessary,  techniques  similar  to  those  described  in 
Malm  et  al.  (2013)  or  Marras  et  al.  (2013b, a),  among  others, 
may  be  valuable  options  to  consider  for  the  solution  of  the 
shallow  water  and  Euler  equations.  Research  in  this  direction 
is  still  on-going;  especially  so  in  atmospheric  applications. 

Finally,  to  assess  the  ability  of  the  algorithm  to  handle 
conforming  and  nonconforming  grids  with  refinement,  we 
complete  the  study  by  building  six  statically  adapted  grids. 
We  constructed  a  hxed  high-resolution  grid  core  in  the 
region  where  the  zonal  jet  is  most  likely  to  develop,  and 
coarsened  the  rest  of  the  domain  with  a  1-  and  2-level  de¬ 
refinement  approach  to  estimate  the  error  when  the  total 
number  of  grid  points  is  drastically  reduced.  The  DG  solution 
algorithm  on  the  refined  grids  is  based  on  the  procedure  by 
Kopera  and  Giraldo  (2013)  and  references  therein. 

The  paper  is  organized  as  follows.  In  Sec.  2,  we  report  on 
the  construction  of  high-order  grids  on  the  sphere.  The  model 
equations  are  described  in  Sec.  3  whereas  the  numerical 
method  of  solution  is  described  in  Sec.  4.  We  describe  the 
tests  and  analyze  the  results  in  Sec.  5.  We  hnally  summarize 
our  hndings  in  Sec.  6. 

2.  High  order  grid  generation  on  the  sphere 

Grid  generation  on  the  sphere  is  a  relatively  easy  task. 
However,  hnding  the  grid  that  can  suit  different  numerical 
methods  is  not.  Because  element-based  Galerkin  methods 
have  the  advantage  of  being  flexible  with  respect  to  the  grid, 
we  analyze  how  different  meshes  on  the  sphere  can  affect  the 
spectral  element  and  discontinuous  Galerkin  solution  of  the 
shallow  water  equations. 

A  standard  approach  to  high  order  grid  generation  is  based 
on  the  construction  of  a  linear  grid  first  (i.e.  a  grid  with 
straight  edges  whose  extrema  are  the  only  points  that  he  on 
the  geometry),  and  then  to  populate  it  with  the  high  order 
internal  points.  The  population  step  is  performed  element¬ 
wise  in  the  logical  space  obtained  by  a  proper  projection 
from  the  physical  space.  Once  the  higher  order  grid  points 
have  been  built  on  the  plane,  a  backward  projection  onto  the 
physical  geometry  moves  the  new  high  order  elements  onto 
the  sphere.  This  process  is  snch  that  the  high  order  elements 
approximate  the  surface  faithfully.  The  projection  from  the 
sphere  to  the  logical  space  may  differ  from  grid  to  grid. 
The  gnomonic  projection  by  Ronchi  et  al.  (1996)  was  used 
to  build  the  high  order  grids.  In  what  follows,  we  describe 
how  the  RLL,  Hex,  and  Ico  grids  were  constructed. 

2.1.  Reduced  Lat-Lon  grid  (RLL) 

Since  the  effort  of  Phillips  (1957)  to  reduce  the  singularity 
problem  at  the  poles  of  a  global  latitude-longitude  grid, 
different  types  of  reduced  Lat-Lon  grid  (RLL  from  now 
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on)  have  been  used.  Partially  based  on  the  composite 
methods  of  Starius  (1977)  and  Browning  et  al.  (1989), 
Lanser  et  al.  (2000)  solved  the  shallow  water  equations 
on  a  Lat-Lon  grid  away  from  the  poles  combined  with  a 
stereographic  grid  at  the  poles.  In  this  paper,  we  build 
and  test  a  modihed,  high  order  version  of  Lanser’s  grids. 
The  interfaces  between  the  reduced  Lat-Lon  area  and  the 
polar  caps  are  obtained  by  a  transhnite  interpolation  (TFI) 
(Gordon  and  Hall  1973;  Eriksson  1984).  The  way  this  is 
done  will  be  described  shortly.  To  build  the  RLL  grid  we 
used  a  multiblock  approach  (see,  e.g.,  Thompson  (1987)) 
that  consists  of  building  a  set  of  independent,  simply 
connected  surface  grids  that  are  eventually  patched  together 
to  form  the  global  mesh.  The  main  RLL  region  and  the 
polar  caps  are  the  first  blocks  to  be  built.  The  main  Lat-Lon 
region  is  composed  of  four  faces  obtained  from  the  master 
face  7i  =  [-7r/4  <  A  <  7r/4]  x  [ifrmn  < ‘P  <  Pmax]  and  from 
its  rotation  with  respect  to  the  z-axis  of  the  sphere.  The 
variables  A  and  ip  indicate  longitude  and  latitude.  The  RLL 
band  in  the  northern  and  southern  hemispheres  is  delimited 
by  Pmin  and  Prnax  ■  The  rotation  matrix  from  71  to  the  three 
remaining  faces  72,73,74  is  obtained  by  the  combined  effect 
of  a  translation  on  the  a:y-plane  and  a  rotation  around  the 
2-axis;  this  transformation  is  given  by 
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The  first  matrix  produces  a  translation  [/  m  n]  =  [—  1  —  1  — 
1]  along  (a:,  y,  z)  and  is  followed  by  the  rotations  \rot,i=2,-i,A  = 
(7r/2,  TT,  37r/2)  given  by  the  action  of  the  second  matrix.  The 
coordinates  on  each  rotated  planar  face  are  simply  a  function 
of  (a:,  y,  2)-^^  and  are  computed  as 

[a;  y  2  1]^,  =  [a;  y  2  l]^i[T],  j  =  2,3,4.  (1) 

Once  rotated,  each  gridded  face  is  projected  onto  the 
sphere.  The  polar  caps  are  built  in  the  same  way,  except 
that  the  starting  point  is  the  square  plane  dehned  within  the 
region  79  =  [— 7r/4  <  A  <  7r/4]  x  [— 7r/4  <  A  <  7r/4].  At  this 
point  we  have  a  linear  grid  made  of  6  disjoint  regular  patches. 
The  top  view  of  this  is  shown  in  Fig.  1.  Eight  new  surfaces 
must  now  be  built  to  hll  the  ungridded  gaps. 


Figure  1.  Linear  grid.  Top  x  —  y  view  of  the  disjoint  lateral  Lat-Lon 
and  top  patches  in  the  multi-block  RLL  grid.  TFI  is  used  to  connect 
them  together  through  an  interface  surface  grid. 


To  build  each  interface  surface  by  TFI  we  need  information 
from  its  four  boundary  edges.  With  reference  to  the 
schematic  of  Fig.  2,  the  region  delimited  by  edges  1,2,  3, 4 
is  built  in  two  steps:  (a)  given  the  coordinates  of  the 

corner  points  1-2-3-4,  build  edges  3  and  4  and  subdivide  them 
into  ID  linear  elements,  where  is  the  user- 

dehned  number  of  elements  along  the  latitude  direction  in  the 
interface  region;  (b )  compute  the  grid  points  in  the  interior  of 
the  patch  using  the  planar  and  linear  transhnite  interpolation 
dehned  by  the  Boolean  sum 

X(|,7))  =U-bV-U®V,  (2) 

where,  given  the  arrays  f  =  0, . . . ,  1  and  r)  =  0, . . . ,  1  in 
computational  space  along  the  local  (u,  v)  directions,  the 
univariate  interpolations  and  tensor  products  (®)  are 
computed  as 

U(|,f))  =  (l-|)X(0,77)  +  |X(l,r)) 
V(|,f))  =  (l-7)X(|,0)+7)X(|,1) 

U  ®  V(|,  77)  =  (1  -  0(1  -  ^)X(0, 0)  +  (1  -  |)77X(0, 1) 

-(1-7))|X(1,0)-|7?X(1,1).  (3) 

In  Equations  (2),  X  is  the  array  of  the  (A,  (^)  coordinates 
of  the  grid  points  in  the  interior  of  the  surface, 
given  the  known  values  of  the  surface  boundary  edges 
stored  in  X(0,  ?)),  X(l,  r)),  X(0  0),  X(0  1),  and  corners 
X(0,0),X(1,1),X(1,0),X(0,1). 


The  high-order  Legendre-Gauss-Lobatto  (LGL)  points  are 
added  to  the  linear  grid  by  projecting  the  linear  elements 
onto  the  auxiliary  gnomonic  space  (on  the  plane),  and  back- 
projecting  the  new  high  order  quads  onto  the  sphere.  The 
total  number  of  elements  and  grid  points  Np  of  the  high 
order  conforming  grid  are  given  by: 
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(5) 


In  (4)  and  (5),  is  the  number  of  elements  of  order  N 

along  the  longitude  direction  of  one  face  in  the  main  Lat- 
Lon  band  that,  by  construction,  coincides  with  the  number 
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of  longitude  elements  in  the  interface  and  polar  regions. 

and  ^  are,  respectively,  the  number  of  latitude 
elements  of  the  Lat-Lon  region  and  of  the  interface  patches. 
The  number  of  latitude  elements  in  the  polar  caps  is  also 
given  by 

Fig.  3  shows  an  example  of  a  high-order  conforming  RLL 
tessellation. 


Figure  3.  Conforming  RLL  grid. 


2.3.  Quad-based  icosahedral  grid 

The  quad-based  icosahedral  grid  of  order  N  is  derived  from 
an  initial  icosahedron  of  20  triangles.  For  better  properties  of 
the  high  order  grid,  every  triangle  is  mapped  onto  a  gnomonic 
space  by  the  transformation  formulas  (Giraldo  (2001)) 

a;  =  rtanA',  (6a) 

?/ =  rtan^s'sec  A',  (6b) 

where  r  is  the  radius  of  the  sphere  and,  given  the  centroids 
(Ac,  ipc)  of  each  triangle,  we  dehne 


A^  =  atan 


cos((9sin(A  —  Ac) 
cos  (pc  sin  (fi  -I-  cos  (fc  cos  tp  cos(A 


(7a) 


p'  =  asin  [cos  pc  sin  p  —  sin  pc  cos  p  cos(A  —  Ac)] .  (7b) 

After  mapping,  the  triangles  are  subdivided  into  smaller 
ones  by  a  Lagrange  polynomial  of  order  n/. 

The  number  of  quadrilateral  elements  and  grid  points  of 
the  final  grid  are  then  given  by 


Np  =  6(iVj  -  2)p2  +  2 

TVc  =  6(AfJ  -  2) 


2.2.  Cubed  sphere  (or  Hexahedral) 

The  cubed  sphere  (see,  e.g.  Ronchi  et  al.  (1996);  Taylor  et  al. 
(1997,  2013))  is  derived  from  the  projection  of  a  cube  onto 
the  sphere.  As  such,  it  consists  of  6  faces  that  are  then 
subdivided  into  as  many  elements  as  necessary.  Like  RLL, 
we  built  this  grid  in  a  multi-block  fashion  by  going  through 
the  same  steps  described  above.  The  fundamental  difference 
is  that  the  hexahedral  grid  (Hex,  from  now  on)  has  only  6 
faces  and  they  are  all  equal.  An  example  is  shown  in  Fig. 
4.  The  internal  LGL  points  are  omitted  for  simplicity.  By 
construction,  the  elements  that  are  furthest  from  the  center 
of  each  face  are  smaller  than  those  at  the  face  center;  in 
addition,  their  distorted  shape  prove  to  be  a  concern  when 
the  resolution  is  not  sufficiently  high.  We  will  discuss  this 
issue  shortly. 


where  the  number  of  points  Np  of  the  original  triangular 
grid  can  be  found  in  Giraldo  (2001). 

Fig.  5  shows  the  8*^  order  icosahedral  grid.  As  for  the 
previous  cases,  the  elements  are  curved  and  lie  on  the 
spherical  surface. 


3.  Shallow  water  equations  on  the  sphere 
3.1.  Equations  of  motion 

Under  the  assumption  of  a  shallow  depth,  the  motion  of 
an  incompressible  fluid  may  be  described  by  the  non-linear 
shallow  water  equations.  The  shallow  water  assumption  holds 
as  long  as  the  characteristic  horizontal  extension  of  the  fluid 
is  much  larger  than  its  depth.  The  vector  representation 
of  the  shallow-water  equations  in  the  Cartesian  coordinates 
X  =  {x,  y,  z)  is  given,  in  flux  form,  by 
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—  +  V  •  {(pu)  =  0  (8a) 

— h  V  •  (0u  ®  u)  =  —(j)V(j)  —  /  (x  X  0u)  +  fix  +  6i^V‘^{(j)u} 

*  (8b) 

where  (p  =  gh  is  the  geopotential  height  (g  and  h  are  the 
modulus  of  the  acceleration  of  gravity  and  the  vertical 
height  of  the  fluid),  vV^  is  the  artihcial  viscosity  term  of 
viscous  coefficient  z/  =  1  x  10®  m^  s“^  that  is  de-activated  by 
the  switch  <5  =  0,  u  =  (zt,  w,  w)^  is  the  velocity  vector,  V  = 
ipx,  dy,  dz)^ ,  and  /  =  2u)zlr^  is  the  Coriolis  parameter  with 
Earth’s  angular  velocity  a;  =  7.292  x  10“®  s“^  and  mean 
radius  r  =  6.37122  x  10®m.  For  the  fluid  particles  to  remain 
on  the  spherical  surface  defined  in  3D  cartesian  space,  the 
normal  component  of  the  velocity  held  with  respect  to 
the  sphere  is  removed  by  the  term  fix  in  the  momentum 
equation,  where  fi  is  the  Lagrange  multiplier  that  serves 
this  purpose.  The  construction  of  the  Lagrange  multiplier 
is  described  in  Sec.  3.2. 

3.2.  Forcing  the  fluid  onto  the  sphere  by  a  Lagrange 
multiplier 

We  used  the  approach  of  Cote  (1988).  Specihcally,  at  the  end 
of  every  time  step  n  +  1,  the  momentum  of  the  how  must  be 
corrected  for  the  particles  to  remain  on  the  sphere.  If  c  and  u 
denote  the  constrained  and  unconstrained  values  of  0u,  the 
new  momentum  is  corrected  according  to 

((^u)”+^  =  (<('U);)+^  +  fix.  (9) 

For  the  huid  to  remain  on  the  sphere,  the  condition 
u  •  X  =  0  must  hold  (i.e.,  the  velocity  held  and  the  position 
vector  must  be  orthogonal),  hence  implying  that  ^Uc  •  x  =  0 
in  equation  (9).  For  this  to  be  true,  the  value  of  the  multiplier 
is  found  to  be  /i  =  — •  x/r^. 

The  Galerkin  discretization  of  (8)  is  described  in  the  next 
section. 

4.  Numerical  method 


where  F(q)  indicates  the  numerical  hux  of  the  advection 
and  diffusion  terms  across  the  boundary  F.  In  the  case  of 
spectral  elements,  where  the  basis  functions  are  continuous 
across  neighboring  elements,  and  given  that  the  problem 
is  solved  on  a  closed  manifold,  the  boundary  integral 
vanishes.  On  the  other  hand,  in  the  case  of  the  DG  method, 
the  boundary  integral  remains  because  of  the  non-zero  inter¬ 
element  fluxes.  Furthermore,  the  surface  integrals  for  DG 
are  dehned  element-wise  over  fig  rather  than  globally  over  fl. 
This  makes  Equations  (11)  and  (12)  a  system  of  equations 
that  are  coupled  via  the  flux  of  Equation  (12).  Both  options 
are  available  within  the  same  code  used  for  this  study  and 
both  methods  share  most  of  the  numerical  machinery. 

The  integrals  above  are  solved  element-wise  on  the  discrete 
polyhedral  approximation  .  The  discrete  domain  is 
dehned  by  the  union  of  N^.  non-overlapping  quadrilateral 
elements  fig.  For  every  element  we  dehne  a  bijective 
transformation  Fqh  :  {x,y,  z)  ^  {^,rf.,Q  that  maps  the 
physical  coordinate  system,  {x,  y,  z)  to  the  reference  system 
{i,g,C)  and  is  snch  that  the  reference  element  Cl’^{^,rf)  = 
[—1,1]  X  [—1,1]  lies  on  the  spherical  surface.  C  is  thus  the 
spherical  radius  vector  whose  discrete  values  identify  a 
spherical  shell  (e.g.  C  =  1  is  the  sphere  of  unit  radius.) 

The  Jacobian  matrix  of  the  transformation  has  com¬ 
ponents  =  d^F’'  and  determinant  j  Jj  =  d^x  ■  G,  where 
G  =  d^x  X  dyX  is  the  surface  conforming  component  of  the 
transformation  (Song  and  Wolf  1999;  Giraldo  2001). 

Given  the  dehnition  of  the  reference  element,  the  element¬ 
wise  solution  can  be  approximated  by  the  expansion 

{N+lf 

9iv(x,t)|ne  =  X!  V’fc(-Fs^^^(x))g7v(xfc,t),  e  =  l,...,fVg 

k=l 

(13) 

where  {N  +  ifl  is  the  number  of  collocation  points  within  the 
element  of  order  N,  and  ipk  are  the  interpolation  polynomials 
evaluated  at  point  k.  The  basis  functions  ipk  are  constrncted 
as  the  tensor  product  of  the  one-dimensional  functions 
hi{^(x))  and  hj{ri(x))  as: 


4.1.  Spectral  elements  and  discontinuous  Galerkin 
approximation 


To  solve  the  shallow  water  equations  by  element-based 
Galerkin  methods  on  a  domain  fl,  we  proceed  by  dehning 
the  weak  form  of  (8)  that  we  hrst  write  in  compact  form  as 


^+V-F(q)  =  S(q),  (10) 

where  q  =  [(p,  <^u]^  is  the  array  of  the  solution  variables,  and 
F  and  S  are  the  flux  and  sonrce  terms  that  can  be  easily 
identified  in  the  eqnations  (8)  above.  Given  a  basis  function 
ip  that  belongs  to  the  Sobolev  space  in  the  case  of  spectral 
elements,  and  to  LS'  in  the  case  of  DG,  the  weak  form  of  (10) 
is  given  by  the  integral 


iP 


g+VF(q)-S(q) 


dLl  =  0, 


(11) 


=  hr{£,{x))  ®  hj{rf{x)), 

\/i,j  =  0, ...,  N.  hi{^{x))  and  hj{rf(x))  are  the  basis  functions 
associated  with  the  N  LGL  points  and  rfj,  respectively, 
given  by  the  roots  of 


(1  -  e)PN{o  =  0, 

where  P'pf(0  is  the  derivative  of  the  At*^-order  Legendre 
polynomial.  Given  these  definitions,  the  one-dimensional 
Lagrange  polynomials  hfl^)  are 


hriO 


1  (1  -  e)p'Ao 

iV(7V+l)(?-C»)Piv(e.)' 


The  functions  hj{g)  are  computed  in  the  same  way. 

The  same  expansion  above  is  used  to  construct  the 
derivatives.  We  write: 


that,  after  integration  by  parts  of  the  flux,  becomes 


/  pj^dn-  /  Vz/>-F(q)dn-  /  ipS{q)dn  + 

Jn  Jn  Jn 

/  z/.n  •  F(q)  dF  =  0,  (12) 


i9gjv(x,t) 

dt 

clgAr(x,t) 
dx  ' 


(N+lf 


clg]v(xfc,t) 


k=l 


=  E 


^V’fc(_£^(x)) 

dx 


dt 


gNi^kfl)-  (14) 
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Mar r as  et  al. 


The  integrals  are  computed  by  numerical  quadrature  on 
the  reference  element  as  follows: 


ini 


q(x,t)dx=  /  q(^,t)|J(0|d^  ■ 


ini 


Nu 


i,j=i 


(15) 


where  the  Gaussian  quadrature  weights  Wij  are  computed  as 
N{N+l)  ■ 

In  the  case  of  spectral  elements,  the  substitution  of  the 
expansions  (13,14)  into  the  weak  form  (12)  yields  the  semi¬ 
discrete  (in  space)  matrix  problem 

^=D^F(q)  +  S(q)  (16) 

where,  for  the  global  mass  and  differentiation  matrices,  M 
and  D,  we  have  that  D  =  M“^D.  The  global  matrices 
are  obtained  from  their  local  counterparts,  and  D®, 
by  means  of  the  direct  stiffness  summation  operation  that 
maps  the  local  degrees  of  freedom  of  an  element  fig  to  the 
corresponding  global  degrees  of  freedom  in  and  adds 
the  element  values  in  the  global  system.  By  construction,  Af 
is  diagonal  (assuming  inexact  integration),  with  an  obvious 
advantage  if  explicit  time  integration  is  used. 

Concerning  the  discontinuous  Galerkin  approximation,  the 
problem  at  hand  is  solved  only  locally,  and  unlike  the  case 
of  spectral  elements,  the  flux  integral  in  Equation  (12)  must 
be  discretized  as  well.  The  element-wise  counterpart  of  the 
matrix  problem  (16)  is  then  written  as: 


^  =  -(M^-^)^F(q'=)  +  (D‘=)^F(q'=)  +  S(q'=),  (17) 

where  we  obtain  M®’®  =  from  the  element 

boundary  matrix,  and  the  element  mass  matrix, 

Af®.  Out  of  various  possible  choices  for  the  dehnition  of 
the  numerical  flux  F(q)  in  Equation  (17),  we  selected  the 
Rusanov  flux  that,  for  the  inviscid  part  of  the  flux,  is 
constructed  as 


where  R(q)  indicates  the  right-hand  sides  of  equations 
(16,17),  K  is  the  number  of  stages  of  the  RK  method,  and 
qO  _  qra  qif  _  qii-Hi  _  computations  were  executed 

with  a  maximum  Courant-Friedrichs-Lewy  number  (CFL, 
Courant  et  al.  (1967))  approximately  0.5.  This  value  ensures 
the  stability  of  the  scheme  given  that  the  time  step  At  is 
such  that 


At  <  CFL  min  [|u  •  x\  +  4>^/X  ■  x] 

xeoh 

The  local  grid  distortion,  is  dehned  as  the  vector 

^  \Ai  At]'  Ai  At]'  Ai  At]  ) 

where  and  Arf  are  the  local  average  grid  size  across  the 
surface. 

In  the  absence  of  dissipation  (e.g.,  artificial  diffusion),  the 
high  order  numerical  solution  of  advective  problems  may  be 
characterized  by  spurious  high  frequency  modes.  To  preserve 
full  stability,  along  with  respecting  the  CFL  condition,  we 
applied  the  standard  spatial  Boyd-Vandeven  filter  (Boyd 
1996)  at  every  time  step.  The  filter  constant  was  chosen  to 
reduce  by  5%  the  highest  modes  only. 

5.  Tests 

To  evaluate  the  behaviour  of  the  solution  on  the  three 
grid  families,  the  highly  nonlinear  zonal  dynamics  test  by 
Galewsky  et  al.  (2004)  was  run  at  different  resolutions  using 
elements  of  order  8.  Based  on  the  analysis  by  Gabersek  et  al. 
(2012),  8  is  the  order  that  we  are  likely  to  use  when  running 
operational  (weather  prediction)  simulations  with  NUMA. 
The  test  consists  of  (i)  an  unperturbed  mid-latitude  zonal 
jet  on  top  of  a  geostrophically  balanced  geopotential  height, 
and  of  (ii)  its  perturbed  counterpart  where  the  perturbation 
is  triggered  by  a  bump  in  the  geopotential.  In  either  case,  the 
initial  zonal  velocity  and  height  are  a  function  of  the  latitude 
ip  (measured  in  radians  unless  otherwise  specified)  as 


u{ip) 


f « 

1  ^!^^exp 

1 

ii  tp  <  fo  OT  (f  >  ipi 
if  po  <ip  <  pi, 

(19) 


F(q)  =  ^  [F]v(q^)  +  FAr(q^)  -  |A|(qS  -  q^)n]  ,  (18) 

where  |A|  =  |n  •  u| -|- is  the  maximum  wave  speed  of 
the  shallow  water  equations.  By  identifying  an  intrinsic 
tangential  direction  in  each  inter-element  boundary  edge,  q^ 
and  q^  denote  the  solution  in  the  left  (L)  and  right  (R) 
elements,  respectively.  For  a  detailed  description  of  the  DG 
algorithm  that  was  used  for  this  study,  the  reader  is  referred 
to  Kopera  and  Giraldo  (2013). 


where  itmax=80ms“^  is  the  maximum  zonal  velocity 
and  {po,pi)  =  (7r/7, 7r/2  —  pq)  are  the  lowest  and  highest 
latitudes  that  delimit  the  jet.  At  the  jet  mid-point,  where 
p  =  7r/4,  the  non-dimensional  parameter  €„  =  exp[— 4/((^i  — 
Po)^]  normalizes  the  jet  magnitude  to  Umax-  From  the  initial 
zonal  flow  given  by  (19),  the  height  g(j)  is  found  by  solving 


gH^p)  =  g<t>o 


f  Fu{p 


tan^ 

r 


dp', 


(20) 


4.2.  Time  integration 

The  solution  of  the  systems  of  ordinary  differential 
equations  (16,17)  is  advanced  in  time  by  a  third-order 
strong  stability  preserving  Runge-Kutta  method  (SSP- 
RK3)  (Cockburn  and  Shu  2001;  Spiteri  and  Ruuth  2002) 
that  yields  the  solution  at  step  n  -|-  1  as 

qfc  =  agq"  -p  af  q*^-!  -p  At  R(q'=“^).  for  fc  =  1,  .  .  .  ,  AT, 


where  the  integral  on  the  right-hand  side  is  computed  by 
numerical  quadrature.  As  in  the  reference,  pq  is  chosen 
to  give  a  global  mean  layer  depth  of  10  km.  The  height 
perturbation  of  test  (ii)  is  the  bump  whose  shape  is  given 
by 

4>'{\,p)  =  ^cos  (v5)exp[-A^/Q^]exp[-(((92  -  pif  / 

for  —  TT  <  A  <  TT,  (21) 
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where  A  is  the  longitude,  (j)  =  120  m,  ip2  =  a  =  1/3 
and  P  =  1/15.  The  definition  of  this  bump  is  such  that  the 
perturbation  is  forced  to  zero  at  the  poles. 

The  flows  in  (i)  and  (ii)  are  simulated  along  a  time 
span  of  144  hours.  In  the  absence  of  any  perturbation,  the 
balanced  field  and  zonal  flow  are  expected  to  maintain  their 
initial  state  indefinitely.  The  analytical  initial  field  is  used 
to  compute  the  L2  error  norms  of  the  numerical  solutions  as 
time  evolves. 

Although  its  solution  seems  trivial,  this  test  is  especially 
demanding  for  a  shallow  water  model  because  of  the  large 
sensitivity  of  the  solution  to  the  grid  resolution  and  geometry. 
Unlike  other  tests  that  may  be  more  forgiving  in  this  respect 
(e.g.  the  steady-state  nonlinear  zonal  geostrophic  flow  by 
Williamson  et  al.  (1992)),  when  the  grid  is  not  sufficiently 
fine,  grid  related  oscillations  may  spoil  the  equilibrium  with 
devastating  consequences.  This  is  clearly  visible  in  Fig.  6, 
where  the  initially  unperturbed  and  balanced  zonal  flow 
completely  loses  its  initial  state  when  the  test  is  executed 
on  the  low  resolution  hexahedral  and  icosahedral  grids.  All 
the  runs  were  executed  at  the  equivalent  resolutions  AA  w 
Aip  =  {0.625°,  1.25°,  2.5°,  5.0°}  and  for  a  varying  number  of 
grid  points  Np  =  {12000,25000,50000,100000}.  The  reason 
for  the  two  sets  of  runs  will  be  clarified  in  Remark  1. 

Below,  the  keywords  Equilibrium  and  Perturbed  jet  are 
used  to  identify,  respectively,  problems  (i)  and  (ii). 

Remark  1:  equivalent  resolutions  A  few  words  are  in 
order  regarding  the  issue  of  equivalent  resolution  when  grids 
of  a  very  different  nature  are  compared  against  each  other. 
In  the  case  of  spectral  elements  that  rely  on  LGL  nodes,  the 
distance  between  two  consecutive  grid  points  changes  within 
the  same  element.  Because  of  this,  the  equivalent  angular 
resolution  along  a  latitude  circle  subdivided  into 
elements  of  order  N  is  2'K/(Ne^xN).  Given  this  definition, 
we  add  that  two  resolutions  are  approximately  equivalent 
when,  within  the  region  of  major  interest  (e.g.,  where  the 
dynamics  is  most  likely  to  develop  from  a  certain  initial 
state),  the  size  of  the  angular  resolution  of  two  different 
grids  are  comparable.  In  other  words,  when  we  say  that 
the  resolution  is,  e.g.,  AA  =  0.625°,  it  means  that  the  mean 
angular  distance  between  two  meridians  in  the  proximity 
of  the  jet  measures  0.625°.  This  definition  makes  the 
measurement  on  the  icosahedral  grid  not  straightforward 
because  of  the  irregular  orientation  of  the  elements.  To 
obviate  this  problem,  the  error  norms  in  the  analysis  below 
are  also  compared  with  respect  to  the  total  number  of  grid 
points. 


5. 1 .  Equilibrium 

Fig.  6  shows  the  vorticity  field  A  after  6  days  for  the 
unperturbed  jet  problem  on  the  three  different  conforming 
grid  types  (rows)  and  Np  =  {25000, 50000, 100000}  grid 
points  (from  left  to  right).  For  two  out  of  the  three  grid 
geometries,  the  effect  of  the  misalignment  is  still  important 
for  as  many  as  25000  grid  points,  which  approximately 
corresponds  to  an  equivalent  resolution  Aip  =  1°  on  the 
cubed-sphere.  The  solution  on  the  RLL  grid  can  maintain 
the  jet  with  approximately  80%  fewer  grid  points  than 
the  hexahedral  grid  and  50%  fewer  than  the  icosahedral 
grid.  This  is  not  reason  enough  to  draw  conclusions  on  the 
superiority  of  the  RLL  grid;  rather,  it  should  only  suggest 
that,  unless  a  sufficiently  high  resolution  is  used  with  grids 


Table  1.  Equilibrium.  Normalized  L2  error  of  {(j>,Us,Vs)  at  day  6 
using  CG.  The  error  is  computed  with  respect  to  the  analytical 
initial  fields  given  by  equations  (19)  and  (20),  and  is  reported  against 
the  total  number  of  grid  points. 


CG,  L2 

N  points 

100000 

50000 

25000 

12000 

RLL 

4.878e-8 

1.802e-7 

2.011e-5 

1.466e-4 

Us 

5.910e-5 

1.586e-4 

5.162e-3 

2.701e-2 

Vs 

5.910e-5 

1.586e-4 

5.167e-3 

2.706e-2 

Hex 

6.099e-5 

4.896e-4 

5.440e-3 

5.460e-3 

Us 

6.679e-3 

5.346e-2 

5.458e-l 

5.059e-l 

Vs 

6.679e-3 

5.346e-2 

5.458e-l 

5.059e-l 

Ico 

2.230e-5 

2.593e-4 

1.057e-2 

1.263e-2 

Us 

2.040e-3 

2.313e-2 

8.035e-l 

9.163e-l 

Vs 

2.040e-3 

2.313e-2 

8.035e-l 

9.163e-l 

Table  2.  Equilibrium.  As  table  1,  but  the  error  is  reported  against 
the  equivalent  resolution  AA. 


CG,  L2 


RLL 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

180000 

49000 

16000 

3000 

1.377e-8 

5.017e-7 

1.106e-4 

2.056e-3 

Us 

2.185e-5 

3.167e-4 

2.553e-2 

4.694e-l 

Vs 

2.185e-5 

3.167e-4 

2.555e-2 

4.694e-l 

Hex 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

124000 

31000 

9600 

1500 

5.109e-5 

4.721e-3 

6.336e-3 

1.038e-2 

Us 

5.601e-3 

4.709e-l 

4.832e-l 

6.918e-l 

Vs 

5.601e-3 

4.709e-l 

4.832e-l 

6.918e-l 

that  are  not  aligned  with  the  characteristic  flow  (which  is 
likely  to  be  the  case  in  realistic  simulations),  we  should  not 
trust  the  solution  with  high  fidelity. 

In  Tables  1  and  2  we  report  the  normalized  L2  error 
norms  of  the  geopotential  height  and  zonal  velocity  against 
the  number  of  grid  points  and,  respectively,  the  equivalent 
resolution  AA.  These  errors  were  computed  using  spectral 
elements  on  conforming  grids;  from  here  on,  we  shall 
refer  to  spectral  elements  as  Continuous  Galerkin  (CG)  to 
distinguish  it  from  the  discontinuous  Galerkin  (DG)  method. 

We  mentioned  above  that  the  structure  of  the  icosahedral 
grid  is  somewhat  irregular  and  it  is  hence  not  straightforward 
to  find  a  direct  relationship  between  the  number  of  grid 
points  and  the  equivalent  resolution.  For  this  reason,  we 
report  the  errors  computed  on  the  Ico  grid  in  Table  3 
against  the  equivalent  resolutions  and  against  the  subdivision 
parameters  [Ah  N5  NA  N3  N2  Al]  that  are  defined  by  the 
user  to  construct  the  grid.  We  plot  the  errors  of  Tables  1-3 
in  Figs.  7  and  8. 

Let  us  focus  on  the  curves  with  continuous  lines  as 
they  trace  the  error  on  the  uniform  conforming  grids.  The 
conclusions  drawn  from  the  vorticity  contours  of  Fig.  6  are 
confirmed  by  these  error  estimates.  Furthermore,  from  Fig. 
8  it  is  interesting  to  see  that  the  L2  error  of  the  Hex  solution 
reaches  its  peak  value  at  AA  1.25°.  This  leaves  very  little 
room  for  resolution  options  when  the  cubed-sphere  is  used 
without  viscosity.  The  results  reported  so  far  were  computed 
with  CG  on  conforming  grids.  However,  we  also  ran  the  same 
tests  using  DG  on  every  grid.  As  expected,  the  errors  are 
practically  identical  to  those  obtained  by  CG;  the  solution 
accuracy  for  CG  and  DG  should  be  identical  because  the 
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Hex  25k 


Hex  50k 


Hex  100k 


Rll  100k 


Ico  100k 


Figure  6.  Equilibrium:  Vorticity  field  (xlO^^  after  6  days.  The  solutions  on  the  hexahedral  (top  row),  on  the  RLL  (middle),  and  on  the 

icosahedral  grids  (bottom  row)  are  plotted  for  25000,  50000,  and  100000  grid  points  (from  left  to  right).  The  color  scale  ranges  between  —12  X  10^^ 
s^^  (dark  blue)  and  16  X  10^®  s^^. 


Table  3.  Equilibrium.  As  table  2,  but  on  the  icosahedron.  Because  of  the  geometrical  structure  of  the  icosahedron,  the  correspondence 
between  the  equivalent  resolution  and  the  number  of  subdivision  indicated  by  N6, . . . ,  A^l  is  not  as  straightforward  as  it  is  for  Hex  and 
RLL. 


Geometry 

Approximate  AA 

N6 

0.625° 

N5 

N4 

1.25° 

N3 

N2 

2.5° 

N1 

5.0° 

(t> 

5.064e-6 

2.299e-5 

2.593e-4 

6.157e-3 

1.261e-2 

1.033e-2 

Ug 

4.381e-4 

2.039e-3 

2.313e-2 

5.169e-l 

9.166e-l 

7.774e-l 

Vs 

4.381e-4 

2.039e-3 

2.313e-2 

5.169e-l 

9.166e-l 

7.774e-l 

solution  is  relatively  smooth  (differences  should  only  emerge 
between  the  two  methods  in  the  presence  of  discontinuities). 
For  comparison,  we  report  the  values  of  the  DG  normalized 
error  norms  in  Table  4.  We  will  show  DG  results  on  non- 
conforming  grids  with  static  refinement  in  Sec.  5.3. 

In  Fig.  9  we  plot  the  time  evolution  of  ||^||l2  during  6  days 
for  100000  point  grids.  After  the  gravity  wave  adjustment 
during  the  first  6  hours  (see  Galewsky  et  al.  (2004)),  the 
regularity  of  the  grid  geometry  still  plays  an  important  role  in 
the  evolution  of  the  error  even  at  a  resolntion  that,  in  global 
mode,  is  typically  considered  fairly  high  in  an  operational 
setting. 

So  far,  the  viscous  term  in  the  equations  has  been  ignored 
(i.e.  (5  =  0).  In  the  quest  of  diminishing  the  negative  effects 
of  the  low  resolution  on  certain  grids,  we  turned  viscosity 


Table  4.  Equilibrium.  As  in  table  2,  but  using  DG. 


DG,  La 

RLL 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

180000 

49000 

16000 

3000 

(p 

1.384e-8 

6.556e-7 

1.203e-4 

2.151e-3 

Us 

2.189e-5 

3.261e-4 

2.636e-2 

4.794e-l 

Vs 

2.189e-5 

3.260e-4 

2.640e-2 

4.794e-l 

Hex 

0.625° 

1.25° 

i5° 

5.0° 

N  points 

124000 

31000 

9600 

1500 

(p 

6.649e-5 

4.948e-3 

6.397e-3 

1.540e-2 

Us 

7.291e-3 

4.933e-l 

4.972e-l 

9.369e-l 

Vs 

7.291e-3 

4.933e-l 

4.972e-l 

9.369e-l 
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Figure  7.  Equilibrium  at  day  6:  Normalized  L2  error  norm  of  (f)  (left)  and  Us  (right)  against  the  number  of  grid  points.  The  error  is  computed 
with  respect  to  the  initial  condition  of  the  balanced  problem  on  every  grid. 


Figure  8.  Equilibrium  at  day  6;  Normalized  L2  error  norm  of  (f)  (left)  and  Us  (right)  against  the  resolution  measured  in  degrees  for  the  RLL  and 
Hex  grids.  The  error  is  computed  with  respect  to  the  initial  condition  of  the  balanced  problem  on  every  grid. 


on  and  recomputed  the  normalized  L2  error  norms  that 
we  report  in  Table  5.  We  compare  these  errors  against  the 
inviscid  estimates  using  12000  and  25000  grid  points  and 
point  out  the  instances  when  artificial  viscosity  was  either 
benehcial  or  harmful  to  the  error  measure.  To  interpret  these 
values,  we  need  to  take  a  parallel  look  at  the  contour  plots 
of  Figures  10  and  11.  In  the  case  of  RLL,  where  the  grid  is 
aligned  with  the  characteristic  flow,  it  was  shown  above  that 
the  resolution  does  not  significantly  affect  the  solution.  The 
opposite,  however,  occurs  in  the  case  of  the  other  two  grids. 
By  adding  viscosity,  what  was  a  bad  inviscid  solution  seems 
to  improve  (see  the  flow  pattern  after  6  days  in  Figs.  10,  11, 
and  the  error  measure  in  Table  5);  in  contrast,  the  solution 


that  was  well  behaved  in  the  inviscid  case  has  deteriorated. 
This  behavior  is  expected  for  two  main  reasons:  (1)  the 
viscous  and  inviscid  equations  are,  mathematically,  not  the 
same  set  of  equations,  and  physically  represent  two  different 
dynamics;  (2)  viscosity  is  so  large  that  it  is  damping  all  of  the 
modes  triggered  by  the  grid  on  the  one  hand  but,  in  the  same 
way,  damps  the  solution  also  where  it  should  not.  A  thorough 
analysis  of  the  diffusion  mechanisms  that  are  currently  used 
in  atmospheric  problems  goes  beyond  the  scope  of  this  paper. 
Nevertheless,  it  should  be  emphasized  that  using  an  arbitrary 
diffusion  operator  may  not  actually  improve  the  solution 
quality  even  though  the  final  solution  may  appear  smooth, 
as  we  have  shown  for  the  RLL  grid.  This  simple  result 
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Figure  9.  Equilibrium:  time  evolution  of  the  normalized  L2  error  norm  of  (left)  and  Us  (right)  relative  to  the  initial  condition  computed  on 
100000  grid  points. 


Table  5.  Equilibrium:  viscous  vs.  inviscid  solution.  Normalized  L2 
error  norms  of  {(f>,Us,Vs)  at  day  6  using  CG  with  artificial  viscosity 
(V)  u  =  1  X  10®  m^  s“^  on  12000  and  25000  grid  points.  For  a  direct 
comparison,  the  inviscid  (I)  errors  computed  on  the  same  grids  are 
also  reported,  value  indicates  that  the  use  of  viscosity  induced  an 
error  increase  in  the  solution.  On  the  contrary,  value  indicates  that 
the  error  norm  decreased  by  using  viscosity. 


CG,  L2 

N  points 

12000-V 

12000-1 

25000-V 

25000-1 

RLL 

0 

1.493e-3 

1.466e-4 

1.481e-3 

2.011e-5 

Us 

1.575e-l 

2.701e-2 

1.548e-l 

5.162e-3 

Vs 

1.575e-l 

2.706e-2 

1.548e-l 

5.167e-3 

Hex 

0 

3.490e-3 

5.460e-3 

1.744e-3 

5.440e-3 

Us 

3.171e-l 

5.059e-l 

1.814e-l 

5.458e-l 

Vs 

3.172e-l 

5.059e-l 

1.815e-l 

5.458e-l 

Ico 

<(> 

7.002e-3 

1.263e-2 

3.284e-3 

1.057e-2 

Us 

5.762e-l 

9.163e-l 

2.947e-l 

8.035e-l 

Vs 

5.762e-l 

9.163e-l 

2.947e-l 

8.035e-l 

stresses  the  importance  of  deriving  a  proper  solution-based 
stabilization  mechanism;  we  are  working  on  this  topic  and 
hope  to  report  our  results  in  a  future  publication. 

5.2.  Barotropic  instability 

Unlike  the  equilibrium  problem,  we  do  not  have  an  analytic 
solution  for  the  perturbed  jet.  Because  of  this,  the  solution 
at  the  highest  resolution  AA  =  0.625°  was  considered 
the  reference  to  compute  the  normalized  errors.  At  this 
resolution,  the  RLL  reference  solution  consists  of  180000 
grid  points.  Based  on  the  analysis  of  the  equilibrium  test 
and  the  high  sensitivity  to  the  resolution  in  the  case  of  the 
hexahedral  and  icosahedral  grids,  we  considered  the  high 
resolution  RLL  run  to  be  the  true  solution  against  which 
the  comparison  should  be  performed.  It  is  a  choice  that 
is  also  supported  by  the  fact  that  the  difference  between 


the  RLL  and  Hex  solutions  first,  between  the  Hex  and  Ico 
solutions  then,  and  finally  between  the  solutions  on  the 
RLL  and  Ico  grids,  are  of  the  same  order  of  magnitude 
(approximately  A  =  0(10“®  s“^)).  This  makes  us  confident 
that  the  comparison  with  respect  to  the  RLL  will  not  be 
biased.  After  6  days,  the  instability  has  fully  developed.  The 
vorticity  computed  at  high  resolution  is  plotted  in  Fig.  12. 

There  is  no  visible  difference  among  the  three  solutions 
with  A(p  w  0.625°.  However,  as  the  resolution  is  decreased, 
its  effect  is  much  more  aggressive  in  the  case  of  the  cubed 
sphere  and  icosahedron  using  12000  and  25000  grid  points 
than  it  is  for  RLL.  In  Figs.  13  and  14,  the  structure  of 
the  jet  in  the  unperturbed  region  remains  regular  on  the 
reduced  Lat-Lon  grid,  with  a  limited  effect  onto  the  fully 
developed  jet.  For  the  other  two  grids  to  resolve  the  jet 
properly,  more  than  25000  grid  points  are  necessary,  as  shown 
in  the  comparative  Fig.  15,  where  the  solutions  show  little  to 
no  difference  starting  at  50k  points  when  viscosity  is  turned 
off. 

In  Fig.  16,  we  plot  the  normalized  error  norms  computed 
against  the  total  number  of  grid  points.  Unlike  for  the 
equilibrium  case,  the  error  measures  of  this  problem  are  very 
similar  for  the  three  different  grids.  This  is  expected  because 
of  the  highly  irregular  structure  of  the  jet  which  is  no  longer 
aligned  with  any  grid.  As  the  resolution  is  coarsened,  the 
error  is  sensibly  higher  than  its  equilibrium  counterpart. 
What  can  make  the  difference  in  the  choice  of  the  grid  in 
this  case,  may  be  as  simple  as  selecting  the  proper  number 
of  grid  points  that  keep  the  error  below  a  certain  threshold. 
At  the  coarsest  resolution,  this  is  most  effectively  achieved 
by  the  RLL  grid. 

We  proceed  by  testing  the  implementation  for  noncon¬ 
forming  grids  with  static  mesh  refinement  in  the  section  that 
follows. 

5.3.  Non- conforming  grids 

To  test  the  DG  solver  on  nonconforming  grids,  we  defined 
a  statically  refined  region  in  the  proximity  of  the  zonal 
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Figure  10.  Equilibrium  12000:  Vorticity  field  at  day  6.  Top  row:  inviscid  solution.  Bottom  row:  artificial  diffusion  with  v  =  1  x  lO^m^  s  ^ .  Solution 
on  the  RLL  (left),  Hex  (center),  and  Ico  (right)  grids. 
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Figure  11.  Equilibrium  25k:  like  Fig.  10,  but  using  25000  grid  points. 


wind.  The  high  resolution  static  rehnement  was  positioned 
within  the  latitude  band  35°  <  (pband  <  55°;  this  region 
covers  approximately  10°  above  and  below  the  central-line  of 
the  initial  jet.  Similar  to  the  uniform  case,  we  dehned  a  set  of 
different  equivalent  grid  resolutions  in  the  banded  region.  We 
used  the  values  in  AX^and  ~  ^‘Pband  =  {0.625°,  1.25°,  2.5°}. 
After  fixing  the  resolution  in  this  region,  for  each  grid  we 
added  a  1-level  and  a  2-level  global  (de-)refinement.  The  total 
number  of  grid  points  in  the  global  mesh  is  thus  reduced.  The 
grid  structure  is  such  that  each  rehned  region  has  half  the 
resolution  of  its  neighbors.  If  two  rehnement  levels  are  used, 
and  AA  =  1.25°  characterizes  the  hnest  portion  of  the  mesh, 
the  resolution  furthest  away  from  it  has  a  AA  =  5.0°.  Based 
on  the  errors  on  the  coarse  conforming  grids,  we  did  not  go 
beyond  this  value  anywhere  in  the  domain.  A  total  of  6  non- 
conforming  grids  were  generated  and  are  shown  in  Fig.  17. 

It  is  well  known  that  the  fundamental  advantage  of 
adaptive  grids  (static  and  dynamic)  is  rehected  in  the  lower 


computational  cost  in  terms  of  hoating  point  operations 
since  the  number  of  degrees  of  freedom  of  the  grid  is  lower 
than  its  uniform  counterpart.  In  Tables  6  and  7  we  report 
the  normalized  L2  error  norms  of  {(j),Us,Vs)  as  a  function 
of  the  number  of  points  given  a  certain  resolution  in  the 
region  where  most  of  the  dynamics  happen.  In  particular,  the 
errors  are  computed  for  one  and  two  rehnement  levels.  With 
reference  to  Figs.  7  and  8,  we  observe  that,  for  every  grid,  the 
number  of  grid  points  necessary  to  obtain  a  certain  error  is 
reduced  by  adopting  one  level  of  rehnement.  We  would  expect 
for  the  same  to  be  true  once  the  grid  is  further  coarsened 
using  two  rehnement  levels,  however,  this  is  not  the  case 
in  this  particular  instance.  Due  to  the  non-linearity  of  the 
problem  and  the  motion  of  the  gravity  waves  across  the  whole 
domain  during  the  hrst  hours  of  the  simulation,  the  ehect  of 
the  poorly  resolved  gravity  waves  in  the  southern  hemisphere 
is  clearly  felt  by  the  jet  resolved  at  a  hner  resolution.  Back  to 
the  6  day  error  evolution  (Fig.  9),  this  is  clearly  observed  at 
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Figure  12.  Perturbed  jet.  vorticity  field  (xlO  ^  s  at  day  6.  The  flow  on  the  RLL  grid  (left)  the  Hex  (middle),  and  Ico  grid  (right)  was  computed 
with  Aip  ~  0.625°. 


Figure  13.  Perturbed  jet  12000'.  vorticity  field  (xlO  ^  s  ^)  at  day  6.  Top  row:  inviscid  solution.  Bottom  row:  with  artificial  diffusion  using 
—  1  X  lO^m^  s^^.  Solution  on  the  RLL  (left),  Hex  (center),  and  Ico  (right)  grids. 


all  times.  This  issue  should  not  deter  us  from  using  more  than 
one  level  of  rehnement  in  general;  it  simply  confirms  the  well 
known  fact  that  the  criterion  for  grid  adaptation  should  be 
very  well  matched  with  the  problem  at  hand  across  the  entire 
domain.  For  a  more  thorough  analysis  of  the  the  construction 
of  DG  for  non-conforming  grids,  please,  refer  to  the  recent 
paper  by  Kopera  and  Giraldo  (2013). 

6.  Conclusions 

We  analyzed  the  behavior  of  two  high  order  Galerkin 
methods  -spectral  elements  (CG)  and  discontinuous 
Galerkin  (DG)-  to  solve  the  shallow  water  equations  on  six 
types  of  spherical  grids.  We  adopted  both  conforming  and 
non-conforming  reduced  Lat-Lon  (RLL),  hexahedral  (Hex), 
and  icosahedral  (Ico)  grids.  Among  others,  one  advantageous 
characteristic  of  element-based  Galerkin  methods  is  that  they 
are  not  constrained  by  the  logical  structure  of  the  mesh.  As 
expected,  the  spectral  element  and  discontinuous  Galerkin 
solutions  showed  almost  equivalent  error  measures  when 
executed  in  conforming  mode.  In  the  case  of  non-conforming 
grids,  we  only  showed  the  DG  results  because  we  have  not 
implemented  non-conforming  CG  on  the  sphere.  However,  we 
expect  the  results  for  non-conforming  CG  to  have  the  same 
behavior  as  the  non-conforming  DG  that  was  shown  in  this 
paper. 

St-Cyr  et  al.  (2008)  compared  a  CG-based  shallow  water 
model  on  the  cubed  sphere  to  a  hnite  volume  solver  on 


Table  6.  Equilibrium:  normalized  L2  error  of  {(p,Us,Vs)  at  day  6 
using  DG  on  the  non-conforming,  banded  grids  of  Fig.  17.  The  error 
is  computed  with  respect  to  the  analytical  initial  fields  given  by  Eqs. 
(19,20). 


DG  I-level,  L2 

RLL-NC 

0.625/1.25° 

1.25/2.5° 

2.5/5.0° 

N  points 

43000 

13000 

5000 

<(> 

2.207e-7 

1.548e-5 

1.337e-3 

Us 

1.647e-4 

3.960e-3 

1.734e-l 

Vs 

1.647e-4 

3.960e-3 

1.730e-l 

Hex-NC 

0.625/1.25° 

1.25/2.5° 

2.5/5.0° 

N  points 

31000 

10000 

1500 

<(> 

1.680e-4 

5.210e-3 

8.513e-3 

Us 

1.844e-2 

5.372e-l 

5.783e-l 

Vs 

1.844e-2 

5.372e-l 

5.783e-l 

Ico-NC 

0.625/1.25° 

1.25/2.5° 

2.5/5.0° 

N  points 

61000 

15000 

4000 

<(> 

8.881e-7 

3.241e-3 

1.348e-2 

Us 

8.294e-5 

2.879e-l 

9.729e-l 

Vs 

8.294e-5 

2.879e-l 

9.729e-l 

a  Lat-Lon  grid.  From  their  study,  they  concluded  that 
spectral  elements  are  unable  to  capture  the  solution  at 
certain  resolutions  on  the  cubed  sphere.  In  this  paper  we 
showed  that  the  inaccurate  CG  solution  on  the  cubed  sphere 
is  not  to  be  blamed  on  the  numerical  method  per  se,  but 
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Figure  14.  Perturbed  jet  25k:  like  Fig.  13,  but  with  25000  grid  points. 


Hex  25k 


Hex  50k 


Hex  100k 


Rll  25k 


Rll  50k 


Rll  100k 


Ico  25k 


Ico  50k 


Ico  100k 


Figure  15.  Perturbed  jet:  vorticity  field  (xlO  ^  s  after  6  days.  Top  row:  Hex.  Middle  row:  RLL.  Bottom  row:  Ico.  Solutions  plotted  for  25000, 
50000,  and  100000  grid  points  (from  left  to  right).  The  color  scale  ranges  between  —12  x  10“^  s“^  (dark  blue)  and  16  x  10“^  s“^. 


rather,  on  the  unfortunate  combination  of  the  resolution  and  We  confirmed  that  a  coarse  resolution  is  responsible  for  the 
mesh  alignment  with  the  dynamics  of  the  problem  at  hand,  breaking  of  the  solution;  however,  by  a  simple  change  of 
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Figure  16.  Perturbed  jet  at  day  6:  Normalized  L2  error  norms  of  <f)  (left)  and  Us  (right)  against  the  number  of  grid  points.  The  error  is  computed 
with  respect  to  the  8*^  order  RLL  solution  at  0.625°,  which  corresponds  to  180000  points. 


Figure  17.  High-order  (curved  elements)  non-conforming  grid  with  1-level  (top  row)  and  2-level  (bottom  row)  refinements.  RLL  (left  column), 
Hex  (middle  column),  and  Ico  (right  column).  The  highest  resolution  in  the  main  region  of  interest  is  Ri  Atp  —  0.625°. 


the  grid  geometry,  we  have  disproved  the  responsibility  of 
CG  alone  by  showing  that  the  instabilities  are  merely  an 
undesired  effect  of  the  grid  configuration. 

To  show  the  flexibility  of  high-order  spectral  elements 
and  discontinuous  Galerkin  with  respect  to  the  mesh,  and 
how  the  solution  may  be  sensitive  to  the  grid  configuration, 
we  used  the  highly  non-linear  zonal  flow  tests  proposed 
by  Galewsky  et  al.  (2004).  Galewsky’s  tests  consist  of 
(i)  an  unperturbed  mid-latitude  zonal  jet  on  top  of  a 


geostrophically  balanced  geopotential  height,  and  (ii)  its 
perturbed  counterpart  where  the  perturbation  is  triggered 
by  a  bump  in  the  geopotential. 

The  first  main  conclusion  from  our  analysis  is  that  both 
the  inviscid  and  viscous  solutions  have  shown  an  important 
sensitivity  to  the  mesh  resolution  and  geometry.  Particularly 
so  are  the  solutions  on  the  hexa-  and  ico-sahedral  grids  at 
low  resolution.  We  observed  the  same  behavior  for  both 
Galewskys’s  tests.  In  the  case  of  the  flow  in  equilibrium, 
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Table  7.  As  Table  6,  but  with  2-level  refinement. 


DG  2-level,  L2 

RLL-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

13000 

5000 

0 

1.265e-5 

9.183e-4 

Us 

3.964e-3 

1.599e-l 

Vs 

3.964e-3 

1.599e-l 

Hex-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

9600 

1500 

0 

5.187e-3 

8.398e-3 

Us 

5.364e-l 

6.348e-l 

Vs 

5.364e-l 

6.348e-l 

Ico-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

15000 

4000 

0 

3.670e-3 

1.298e-2 

Us 

3.249e-l 

9.695e-l 

Vs 

3.249e-l 

9.695e-l 

as  the  resolution  is  decreased  below  an  approximate  value 
of  1.25°  in  latitude  and  longitude,  the  misalignment  of  the 
mesh  in  the  Hex  and  Ico  configurations  is  fatal.  The  fact  that 
the  grid  may  not  be  aligned  with  the  flow  (which  is  likely  the 
case  in  most  real-life  simulations)  showed  not  to  be  the  only 
factor  of  concern  in  certain  configurations.  In  fact,  although 
the  flow  is  not  aligned  with  any  of  the  three  grids  in  the  case 
of  the  barotropically  unstable  jet,  the  solution  on  the  reduced 
Lat-Lon  was  only  slightly  affected,  whereas  the  accuracy  was 
completely  lost  for  the  other  two  configurations.  We  execnted 
the  same  tests  with  and  without  artificial  viscosity.  Based 
on  the  results  obtained  with  the  use  of  artificial  viscosity, 
we  believe  that  the  extremely  large  and  unphysical  diffusion 
that  was  added  by  means  of  a  second-order  Laplace  operator 
with  coefficient  i/  =  1  x  10^  m^s“^,  should  be  avoided  for  the 
sake  of  accuracy. 
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